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A 2-D  Boussinesq  model  describing  heat-driven  buoyant  convection  in  a polygonal  en- 
closure is  presented.  The  hydrodynamics  is  based  on  the  time-dependent  Navier-Stokes 
equations  with  constant  viscosity  and  thermal  conductivity;  no  turbulence  model  or 
other  empirical  parameters  are  introduced.  The  polygonal  domain  is  mapped  via  a 
Schwarz- Christoff  el  transformation  onto  a rectangle.  A finite  difference  scheme  second- 
order  in  space  and  first-order  in  time  is  used  to  integrate  the  evolution  equations,  and 
an  elliptic  solver  is  used  to  solve  the  pressure  equation.  Computational  results  for  high 
Reynolds  numbers  are  presented  through  the  use  of  Lagrangian  particles  which  allow 
one  to  visualize  the  flow  patterns. 


1 Introduction 

In  computational  aerodynamics,  calculations  of  the  steady-state  flow  over  aircraft  have  rou- 
tinely required  the  use  of  sophisticated  grid-generating  techniques  [1].  In  other  areas,  such  as 
ship  hydrodynamics,  flow  in  machinery,  and  shore  effects  in  oceanography,  non-rectangular 
grids  are  routinely  used.  In  fire  research,  however,  the  study  of  fluid  flow  in  more  compli- 
cated domains  is  a more  recent  activity.  Many  of  the  techniques  developed  by  the  CFD 
community  for  solving  the  flow  equations  on  non-rectangular  grids  have  one  or  more  of  the 
following  limitations:  1)  They  are  designed  to  compute  only  the  steady-state  solution.  2) 
They  rely  on  various  types  of  turbulence  models.  3)  They  are  inefficient  due  to  the  use  of 
unstructured  grids. 

We  have  previously  developed  an  algorithm  for  the  computation  of  the  time-dependent 
buoyant  convection  induced  by  a room  fire  [2]  - [7].  Although  this  model  has  proven  to  be  a 
useful  tool  in  fire  research,  it  is  limited  to  flows  in  rectangular  enclosures.  In  this  document 
we  generalize  our  algorithm  to  allow  for  polygonally-shaped  domains,  making  it  possible  to 
study  the  geometrical  effects  of  building  elements  such  as  windows,  soffits  and  stairways.  A 
conformal  map  from  the  polygon  to  a rectangle  is  used  to  transform  the  spacial  coordinates 
of  the  Boussinesq  equations.  The  advantage  of  this  approach  is  that  we  retain  the  efficiency 
and  resolution  of  the  original  algorithm.  (Since  no  turbulence  model  is  included,  we  must 
sufficiently  resolve  the  boundary  layer  for  a flow  whose  Reynolds  number  approaches  105,  a 
reasonable  value  for  enclosure  fires.)  The  disadvantage  of  the  conformal  transformation  is 
that  it  can  lead  to  severe  variations  in  grid  cell  size,  and  thus  significantly  limit  the  time-step 
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of  our  explicit  scheme.  For  this  reason,  we  consider  domains  which  do  not  yield  too  large 
a variation  in  the  Jacobian.  We  have  chosen  for  the  moment  to  retain  the  machinery  built 
up  from  our  previous  work  and  implement  the  conformal  transformation  without  changing 
the  basic  algorithm.  However,  for  domains  which  exhibit  a larger  variation  in  the  Jacobian 
(more  than  2 or  3 orders  of  magnitude),  a better  approach  would  be  to  use  an  implicit, 
unconditionally  stable  time-stepping  scheme  [8]. 

The  model  of  buoyant  convection  is  presented  in  Sections  2 and  3 and  a brief  description 
of  the  numerical  method  in  Section  4.  Section  5 contains  a detailed  description  of  the 
algorithm.  Section  6 presents  the  results  of  a few  computations,  including  a verification  of 
the  “trench  effect”  [2]  in  a more  realistic  stairway  geometry.  An  Appendix  is  also  included 
which  contains  the  algorithm  used  to  introduce  and  track  the  Lagrangian  particles,  and  also 
a linear  stability  analysis  of  the  lagged-dissipation  scheme. 


2 Hydrodynamic  Model 

We  consider  a Boussinesq  fluid  with  constant  viscosity  and  thermal  conductivity  in  a polyg- 
onal enclosure  driven  by  a prescribed  heat  source.  We  start  with  the  equations  of  motion 
for  a thermally  expandable  ideal  gas  [3]  in  which  we  include  viscous  dissipation  and  thermal 
conductivity. 
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Here,  all  symbols  have  their  usual  fluid  dynamical  meaning:  p is  density,  ut  are  the  velocity 
components,  p is  pressure,  g is  the  acceleration  of  gravity,  are  the  components  of  the  vector 
describing  the  direction  of  the  gravity  vector,  v is  kinematic  viscosity,  Cv  is  the  constant- 
pressure  specific  heat,  T is  temperature  and  K is  the  thermal  conductivity,  t is  time  and  Q 
is  the  spatially  and  temporally  prescribed  heat  source. 

If  we  combine  these  equations  as  described  in  [4],  nondimensionalize  them  as  described 
in  the  appendix  of  this  reference  and  make  the  Boussinesq  approximation,  we  obtain  the 
following  set  of  equations, 


dul  dui  dp 
dt  + J dxj  ^ dxi 


dp  dp_ 
dt  1 dxi 


- — -Q  + /cV2p 

7 


{p  ~ po(z))ki  = vV2Ui 


duj 

dxi 


= 0 


(2) 


2 


Here,  all  symbols  have  the  meanings  given  above  but  in  dimensionless  form,  po(z)  is  the 
initially  stratified  density  profile  assumed  to  depend  only  upon  the  vertical  coordinate  2, 
and  k is  the  dimensionless  thermal  conductivity.  The  dimensionless  quantities  are  defined 
as  follows:  lengths  are  relative  to  the  height  of  the  enclosure,  pressure  is  relative  to  ambient 
pressure,  time  is  relative  to  the  height  of  the  enclosure  and  the  velocity  scale,  and  velocity 
is  relative  to  a scale  U,  and  the  density  perturbation  is  relative  to  ambient  density  with  a 
small  parameter  /?.  If  we  denote  temporarily  the  dimensional  quantities  by  an  asterisk,  then 
these  relationships  can  be  written  as  follows: 

x*  = Hxi , p*  = po(0)U2p 

t*  = Ht/U , p*  = PO{0){1+ f3p) 

u*  = Uui,  p*(z*)  = po(0)(l  + (3p0{z)) 

U and  (3  are  defined  in  terms  of  the  magnitude  of  the  heat  source  as  follows: 


/?  = U2/(gH) 

U = (Qog/ipoC.ToH^/^iqogKpoC^))1/3 


Here  Q0  is  the  strength  of  the  three-dimensional  heat  source  in  units  of  energy  per  unit 
time,  qo  is  the  strength  of  the  two-dimensional  heat  source  in  energy  per  unit  length  per  unit 
time,  H is  the  height  of  the  enclosure  and  all  other  quantities  have  their  usual  meanings. 
See  [4]  and  [7]  for  more  details  of  the  scaling.  Henceforth,  all  variables  will  be  regarded  as 
dimensionless. 

In  the  two-dimensional  case,  these  equations  can  be  rewritten  as  follows: 
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Note  that  gravity  need  not  be  restricted  to  the  2-direction  (downward),  but  may  point  in 
any  direction  according  to  its  components  gx  and  gz.  The  norm  of  the  gravity  vector  is  unity. 
A no-flux  condition  is  imposed  at  the  boundary.  In  addition,  one  may  choose  a no-slip  or 
a free-slip  condition,  as  well  as  an  adiabatic  or  constant  temperature  condition  at  the  wall. 
The  fluid  is  initially  quiescent. 
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3 The  Transformation 


The  physical  domain  may  be  taken  as  any  simply  connected  polygon  whose  vertices  may 
or  may  not  extend  to  infinity.  This  polygon  is  mapped  conformally  into  a rectangle  whose 
aspect  ratio  is  dependent  on  the  shape  of  the  polygon.  Details  of  the  Schwarz-Christoffel 
transformation  may  be  found  in  [9]  and  a sample  mesh  is  presented  in  Fig.  1.  The  equations 
above  are  written  in  terms  of  the  physical  coordinates  x and  z.  We  shall  take  the  coordinates 
of  the  computational  rectangle  to  be  £ and  p.  Eqs.  (3)  may  then  be  transformed  as  follows: 
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In  transforming  the  viscous  terms  of  the  momentum  equations,  we  have  made  use  of  the 
vector  identity 


V2u  = V(V  • u)  — V x (V  x u) 

and  the  fact  that  V • u = 0.  Because  the  coordinate  transformation  is  conformal,  the 
boundary  conditions  for  Eqs.  (3)  can  be  implemented  in  the  computational  domain  in  exactly 
the  same  way  in  which  they  would  be  for  a rectangular  enclosure. 
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4 Numerical  Methods 


Eqs.  (4)  are  a mixed  parabolic/elliptic  system  of  partial  differential  equations;  be.,  the 
equations  for  the  density  and  for  the  velocity  components  are  parabolic,  whereas  that  for  the 
pressure  is  elliptic.  The  incompressible  equations  of  hydrodynamics  are  well  known  to  have 
this  mixed  character.  Analytical  studies  of  the  ability  of  several  candidate  finite  difference 
schemes  to  calculate  internal  gravity  waves  without  dissipation  [5]  led  to  the  conclusion  that 
methods  of  second  order  accuracy  in  space  and  time  are  highly  desirable  for  systems  of  this 
type. 

We  have  chosen  a lagged-dissipation  scheme  — all  time  derivatives  are  replaced  by  central 
differences  over  twice  the  time-step  size  (a  leap-frog  scheme)  while  all  viscous  and  conduction 
terms  are  differenced  at  the  lagged  time  level.  Other  terms  in  the  evolution  equations  (the 
first  three  of  Eqs.  (4))  are  evaluated  at  the  mid  level  of  the  three  level  scheme.  This  scheme  is 
not  formally  second-order  accurate  in  time.  However,  for  the  Reynolds  and  Prandtl  numbers 
capable  of  capturing  large  scale  features  in  typical  enclosure  fire  simulations  (Re  > 104),  the 
dissipation  terms  are  very  small,  and  the  formal  first-order  accuracy  due  to  dissipation  does 
not  pose  a practical  limitation.  The  grid  is  taken  to  be  uniform  in  the  £ and  rj  directions, 
although  the  spacings  6£  and  Srj  may  be  different.  Within  each  rectangular  cell,  vector 
components  are  evaluated  at  the  sides,  scalar  quantities  at  the  center. 

The  density  evolution  equation  in  continuous  form  is  the  mass  conservation  equation 
minus  the  expression  for  the  velocity  divergence.  Each  of  these  two  equations  is  approximated 
by  central  differences  and  then  subtracted.  The  density  at  all  faces  is  approximated  by  the 
mean  of  the  density  at  the  centers  of  adjacent  cells.  This  procedure  ensures  global  mass 
conservation  as  well  as  second  order  accuracy. 

The  momentum  equations  are  differenced  in  the  vector  invariant  form  shown  in  Eqs. 
(4).  This  ensures  nonlinear  stability  and  complete  compatibility  between  the  “primitive 
variable”  formulation  presented  here  and  a vorticity,  stream-function  formulation  (in  the 
dissipation-free  case  when  the  Jacobian  is  unity),  see  [6]  for  details. 

The  pressure  equation  is  derived  by  differencing  in  time  the  centrally-differenced  in- 
compressibility equation,  and  then  making  use  of  the  discretized  momentum  equations  to 
eliminate  the  time  differences.  The  result  is  an  elliptic  partial  differential  equation  for  the 
pressure.  For  the  Boussinesq  model,  the  linear  algebraic  system  arising  from  its  discretization 
has  constant  coefficients  and  can  be  solved  by  a fast  direct  method,  see  [7]  for  details.  The 
solution  to  the  pressure  equation  constitutes  the  bulk  of  the  numerical  computation  since 
the  density  and  the  velocity  are  updated  explicitly  once  the  pressure  gradients  are  known. 

Stability  of  the  computational  scheme  imposes  a limit  on  the  time-step  size  relative  to 
the  spatial  mesh  size,  see  [4],  [5]  and  Appendix  B.  This  may  be  extremely  limiting  due  to 
the  grid  distortion  brought  about  by  the  conformal  mapping,  whose  Jacobian  may  vary  up 
to  two  orders  of  magnitude.  For  geometries  which  are  even  more  distorting  than  those  we 
present  here,  another  method  of  grid  generation  would  probably  have  to  be  implemented, 
which  would  unfortunately  reduce  the  efficiency  of  the  present  algorithm. 
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5 The  Algorithm 

5.1  The  Density  Equation 

We  now  write  out  the  details  of  the  algorithm.  For  the  density  equation,  return  to  the  full 
continuity  equation  (in  dimensional  form): 


jdp  , dpuj 

dt  d£j 

Use  the  nondimensionalization  presented  earlier 


= 0 


rdp  , A d{p0  + p)  , ri  , of  , ^duj  n 

Jst + [ 0{po  p)]^  ~ 

where  p = 1 + f3(p0  + /?),  where  all  variables  are  dimensionless,  and  where  we  have  divided 
through  by  (3U / H.  If  we  now  formally  allow  /?  — > 0,  then  we  recover  the  Boussinesq  equation 
for  the  density: 

dp  , A d(p0  + p)  , duj  _ 0 


J yrr  + Uj 

dt  d£j 
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The  same  procedure  used  above  for  the  partial  differential  equations  (PDE)  can  now  be 
applied  to  the  finite  difference  equations  (FDE).  We  desire  to  keep  conservation  form  for  the 
FDE  as  well  as  for  the  PDE.  Hence, 
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or,  rewriting, 
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Note  that  the  Jacobian  subscript  p indicates  that  J is  to  be  evaluated  where  p is,  namely  at 
the  center  of  a grid  cell.  This  convention  will  be  followed  from  here  on. 

Now  we  subtract  from  the  density  equation  the  expression  for  the  velocity  divergence 
and  obtain 
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■T  pi,ik  — \Pi+l,k^ik  Pi—l,k'^'i—l,k) 

Fpr),ik  — 08 Tj  Pi,k—l^i,k—l) 

Finally,  use  the  same  nondimensionalization  as  cited  above,  using  pik  = 1 + (3{po,ik  + pik)- 
Note  that  only  p^  depends  upon  time.  Also,  note  that  we  can  divide  by  (3  and  then  allow 
it  formally  to  vanish  as  done  above  for  the  PDE  to  obtain  the  finite  difference  equation  for 
the  Boussinesq  model. 
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Here  k — K/(pCp).  Q will  be  described  below. 

At  boundaries,  the  density  fluxes  are  determined  by  the  no-inflow,  no-outflow  condi- 
tions. Also,  we  must  specify  adiabatic  or  cold-wall  boundary  conditions,  which  determine 
the  temperature  and  hence  the  density  (since  the  perturbation  density  is  the  negative  of  the 
perturbation  temperature)  in  ficticious  cells  adjacent  to  the  boundaries. 
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Here,  the  plus  sign  corresponds  to  adiabatic  boundary  conditions  (BC)  and  the  minus  sign  to 
cold  wall  BC.  Adiabatic  BC  imply  zero  derivative  in  the  temperature  (density)  perturbation 
across  the  boundary  whereas  cold  wall  BC  imply  the  temperature  (density)  perturbation  is 
zero  at  the  boundary. 

The  discretization  of  the  time  derivative  has  been  chosen  as  a leap  frog.  For  the  second- 
order  time  steps 

pT  = Pa1  + 2 « 

and  for  the  first-order  time  steps 


pT  = Pl  + StRn 


p,ik 


IX, 


ik 


where 

Ki * = - K,#  - + K ((vVoU  + (VV)3TJ) 

Rlik  = ~ ~ ^Ql  + * ((vV»U  + (V2p)l) 

The  Laplacian  operator  V2  has  been  defined  above.  Notice  that  we  use  a lagged-diffusion 
temporal  discretization  for  the  dissipation  on  the  second-order  time  steps  only. 

It  is  more  convenient  to  define  the  heat  source  in  terms  of  the  computational  coordinates 
£ and  77,  even  though  it  would  be  preferable  to  use  the  actual  physical  coordinates  x and  z. 
The  trouble  with  the  latter  approach  is  that  it  is  difficult  to  interpolate  a particle’s  position 
in  the  physical  plane,  because  the  mesh  is  obviously  not  uniform  there.  Therefore  in  these 
calculations,  the  following  forms  for  the  heat  source  have  been  used: 

Qik  — fn  Jp,ik  Qt,iQv,k 


where  fn  is  either 
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The  latter  function  represents  a buoyant  thermal.  The  spatial  factors  are 
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5.2  The  Momentum  Equations 

For  the  velocity  equations,  we  use  the  following  difference  scheme: 
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The  derivatives  dxi/d^m  are  evaluated  at  the  corners,  thus  g^ik  is  evaluated  wherever  Uik  is, 
and  gvjk  is  evaluated  wherever  Wik  is.  To  improve  the  computational  efficiency,  we  set 

.ri  &u>,ib  + 1Vu/,i,k—l  ^i,k—l  Pi+l,k  d-  Pik  ^ik  ^i,k— 1 

G(,,k  = j 9iA 1 + i/ Tr, 

x-,  Vu ifik  VJik  d”  i,k  LVi— l,k  Pi,k+ 1 d-  Pik  ^ik  l,fc 

^Vpk  = 2 9r\,ik  ^ ^ ££ 

This  notation  will  be  of  use  in  the  pressure  equation. 

In  order  to  compute  these  terms  at  the  boundary  we  must  either  apply  free-slip  or  no-slip 
boundary  conditions.  For  0 < k < K, 


and  for  0 < i < I 


where  the  plus  sign  corresponds  to  free-slip  and  the  minus  sign  to  no-slip  BC. 

The  temporal  derivatives  are  handled  by  a leap-frog  discretization.  For  the  second-order 
time  steps  using  the  lagged-dissipation  scheme, 
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and  for  the  first-order  time  steps, 
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Here  again,  G and  G are  the  same  except  that  the  former  has  a lagged  dissipation  term, 
while  the  latter  does  not. 


5.3  The  Pressure  Equation 

The  incompressibility  condition  is 


J. 


p,ik 


St 


Sr] 


= 0 


and  this  condition  can  be  used  to  derive  the  linear  algebraic  equation  system  for  the  pressure. 
Away  from  boundaries,  that  is  for  2 < i < I — 1,2  < k < K — 1,  we  have 


Pi  + l,k  2pik  Pi— l,k  Pi,k+1 


2 pik  T Pi,k— 1 G^^ik  l,fe  G^^ik  G-q^^k—  1 


S? 


Srp 


Sp 


(5) 


The  lack  of  a superscript  implies  that  all  quantities  are  to  be  evaluated  at  time  level  n. 
For  cells  adjacent  to  the  left  boundary,  we  have  iiQk  = 0 for  all  n,  thus 


Plk  P0k 

~~sf~  ” 

Subtracting  this  equation  from  the  general  equation  for  pressure  with  i — 1 yields  the 
equation 

p2,k  ~ Pl,k  Pl,k+1  ~ + p\,k-l  _ Q,lfc  GVtik  — GVti,k-l 

St2  Sp2  SS,  Sr] 

In  a similar  fashion,  the  above  equation  can  be  derived  for  the  three  other  edges  and  the 
four  corners.  To  improve  the  computational  efficiency  of  the  code,  we  set  G^,ofc  = G^jk  = 0 
for  1 < k < K,  and  we  set  GVjl o = Gv,iK  — 0 for  1 < i < I.  Then  the  index  range  of  Eq. 
5 may  be  extended  to  the  boundary  cells,  and  we  solve  the  equivalent  Poisson  equation  in 
which  the  normal  derivative  of  the  pressure  is  specified  to  be  zero  at  the  boundary. 
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6 Results 


In  this  section  samples  of  results  generated  using  a code  implementing  the  algorithm  de- 
scribed above  are  presented  and  discussed.  The  results  represent  some  problems  of  interest 
to  scientists  concerned  with  enclosure  fires.  The  computations  were  performed  on  the  IBM 
RISC  System/6000  Model  550  in  the  Computational  Combustion  Facility,  which  is  a joint 
facility  of  the  Computing  and  Applied  Mathematics  Laboratory  and  Building  and  Fire  Re- 
search Laboratory  at  NIST.  The  graphics  were  generated  on  a Silicon  Graphics  Personal  Iris 
4D20,  which  is  also  part  of  this  facility. 

The  scale  of  these  computations  is  substantial.  Runs  have  now  been  performed  using 
more  than  a quarter  of  a million  grid  points;  they  require  at  least  50  megabytes  of  memory 
and  take  up  to  24  hours  of  CPU  time.  They  are  performed  in  64  bit  arithmetic.  Typically, 
the  CPU  time  per  node  per  time-step  is  on  the  order  of  2.0  x 10-5  seconds.  If  we  were 
to  save  any  major  fraction  of  the  data  generated,  we  would  be  overwhelmed.  Therefore,  a 
significant  effort  has  been  expended  in  trying  to  select  only  those  data  required  to  understand 
the  phenomena  being  studied.  We  have  concluded  that  Lagrangian  particle  tracking  for 
transient  phenomena  is  the  most  convenient  method  for  saving  and  examining  data  from  our 
computations,  and  this  procedure  is  described  in  [2], 

The  resolution  of  the  computation  determines  the  maximum  Reynolds  number  of  the 
flow.  Since  the  size  of  the  Reynolds  number  for  the  2-D  flow  scales  as  K 2,  flows  with 
Reynolds  numbers  of  nearly  one  million  can  in  principle  be  computed.  We  have  not  yet  com- 
puted flows  with  this  large  a Reynolds  number,  primarily  because  we  have  been  examining 
non-rectangular,  polygonal  geometries.  For  transformations  such  as  the  stairway,  we  have 
computed  results  for  flows  with  Reynolds  numbers  of  up  to  105,  which  is  reasonably  close 
to  full  scale  for  many  fire  scenarios.  A big  constraint  in  these  computations  is  the  time-step 
size,  which  is  constrained  by  the  CFL  condition  of  the  cells  whose  areas  are  extremely  small. 
Fig.  1 displays  a crude  version  of  the  mesh  used  for  a simulation  of  smoke  through  a window. 
The  Jacobian,  which  is  a measure  of  the  area  of  an  individual  cell,  varies  from  about  10-2 
to  101  for  this  geometry.  Because  the  time-step  size  is  proportional  to  the  square  root  of  the 
Jacobian  (see  Appendix  B),  the  calculation  is  slowed  considerably  by  the  small  cells.  For 
this  reason,  the  present  methodology  is  limited  to  enclosures  whose  cell  areas  vary  no  more 
than  those  shown  in  Fig.  1.  An  alternative  approach  for  difficult  geometries  would  be  to  use 
an  implicit  time-stepping  scheme. 

6.1  The  “Trench  Effect” 

Fires  in  buildings  involve  the  transport  of  heat  and  mass  by  gravity-induced  or  buoyant 
convection.  Generally,  this  convection  occurs  in  rectangular  enclosures  where  the  direction 
of  gravity  is  parallel  to  the  surfaces  of  the  enclosure,  the  walls.  However,  under  certain 
circumstances,  such  as  a fire  in  a stairwell  or  an  escalator,  the  enclosure  may  be  sloped  relative 
to  gravity.  A very  important  example  of  a fire  in  a sloped  enclosure  was  the  devastating 
fire  in  the  King’s  Cross  underground  station  in  England  in  1987,  where  there  was  significant 
loss  of  life  as  well  as  property  damage.  Numerical  simulation  of  this  fire  uncovered  an 
unexpected  phenomenon  which  caused  a very  rapid  spread  of  the  fire  and  led  to  much  of  the 
devastation  [10],  [11].  This  phenomenon  was  termed  the  “trench  effect”,  and  caused  some 


11 


controversy  during  investigations  of  the  King’s  Cross  fire  in  England.  The  phenomenon  was 
ultimately  confirmed  by  experiments  and  additional  simulation,  but  transient  aspects  of  the 
fire  simulation  are  still  of  interest. 

We  repeat  here  an  experiment  performed  with  a rectangular  domain  tilted  35  degrees 
from  the  horizontal  which  was  intended  to  model  a stairwell  [12].  The  major  difference  in 
the  present  calculation  is  the  addition  of  more  realistic  “landings”  to  the  geometry.  Fig.  2 
presents  a time  sequence  for  the  flow  generated  in  an  enclosure  by  a heat  source  (fire)  located 
near  the  base  of  the  stairwell.  The  mesh  size  is  1024  x 256  (262,144  grid  cells),  the  Reynolds 
number  based  on  the  height  of  the  landings  is  4.0  x 104,  and  adiabatic,  free-slip  boundary 
conditions  are  imposed.  The  average  time-step  size  is  8.6  x 10-4  dimensionless  time  units. 
The  CPU  time  per  node  per  time-step  is  2.2  x 10-5  seconds.  In  this  figure,  the  plume  rises, 
but  is  bent  back  toward  the  lower  landing.  After  the  hot  gases  hit  the  ceiling,  they  progress 
both  toward  the  back  wall  and  up  the  ceiling  toward  the  high  end.  However,  the  hot  gases 
leaving  the  heat  source  are  pinned  along  the  floor  and  form  a hot  gas  jet  which  progresses 
up  along  the  floor,  shedding  hot  gases  near  its  front.  This  phenomenon  we  interpret  as  the 
“trench  effect”. 

6.2  Smoke  through  a window 

In  Fig.  3,  we  present  the  flow  induced  by  a fire  in  a small  enclosure  with  a window.  The  mesh 
size  is  256  x 128  (32,768  grid  cells),  the  Reynolds  number  based  on  the  height  of  the  small 
room  is  104,  and  adiabatic,  no-slip  boundary  conditions  are  imposed.  This  geometry  pushes 
the  present  capabilities  to  the  limit.  The  average  time-step  size  for  this  run  is  7.7  x 10-4. 
This  is  due  to  the  extremely  small  cells  around  the  window  soffits.  The  CPU  time  per  node 
per  time-step  is  1.6  x 10-5  seconds. 

As  before  we  notice  that  the  vorticity  of  the  escaping  plume  tends  to  pin  the  hot  gas  to 
the  wall  of  the  enclosure,  a phenomenon  which  has  been  seen  in  actual  fires  where  the  flames 
and  hot  gases  pour  out  of  an  enclosure  opening  and  then  spread  rapidly  up  the  side  of  the 
enclosure.  In  Fig.  4,  we  freeze  in  time  four  calculations  with  different  boundary  conditions. 
There  is  a dramatic  difference  between  the  free-slip  and  the  no-slip  cases.  Ultimately,  all 
demonstrate  the  above  described  phenomena,  however,  the  free-slip  cases  allow  for  a faster 
escape  of  gas  from  the  small  room.  This  type  of  behavior  has  been  seen  in  the  case  of  a 
stairwell  [12],  where  the  “trench  effect”  of  Fig.  2 is  not  as  dramatic  when  no-slip  boundary 
conditions  are  employed. 

6.3  Corridor  with  soffits 

Finally,  in  Fig.  5 we  present  the  movement  of  smoke  down  a corridor  with  uniformly  spaced 
soffits.  The  mesh  size  is  512  x 128,  the  Reynolds  number  is  1.0  x 104,  and  adiabatic,  no-slip 
boundary  conditions  are  imposed.  This  example  points  out  two  difficulties  we  have  expe- 
rienced with  the  method  described  above.  First,  the  computation  was  slowed  considerably 
by  the  small  grid  cells  clustered  around  the  soffits.  This  particular  run  required  about  18 
hours  of  CPU  time;  a similar  run  without  the  soffits  required  about  a tenth  of  the  time. 
Second,  we  observed  that  the  calculation  became  unstable  near  the  end  of  the  run.  Fig.  6 
displays  the  perturbation  density  of  the  boundary  cells,  with  the  enclosure  boundary  serving 
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as  the  spacial  axis.  Where  the  curve  falls  inside  the  enclosure,  the  perturbation  density  is 
negative  and  thus  the  temperature  is  greater  than  its  original  value.  Where  the  curve  falls 
outside  the  enclosure,  the  temperature  is  less  than  its  original  value.  Obviously,  the  density 
profile  at  time  t = 18  reflects  the  instability.  We  suspect  that  the  time-stepping  scheme 
is  responsible,  although  the  linear  stability  analysis  fails  to  identify  the  problem.  Initially, 
it  appeared  as  if  the  boundary  layer  was  not  being  fully  resolved,  but  we  have  noticed  the 
instability  for  Reynolds  numbers  well  below  the  theoretical  limits.  Also,  the  instability  did 
not  develop  until  tens  of  thousands  of  time  steps  had  been  completed,  leading  us  to  suspect 
some  non-linear  effect  of  the  time-stepping  scheme.  We  applied  a second  order  Runge-Kutta 
scheme  to  replace  the  lagged-dissipation  scheme,  and  the  instability  did  not  reoccur.  The 
Runge-Kutta  scheme,  however,  requires  two  calls  to  the  pressure  solver  per  time  step  instead 
of  one,  which  nearly  doubles  the  computing  time.  This  is  the  price  required  to  remove  the 
long-time  instability. 

A Appendix:  Particle  Tracking  Algorithm 

A.l  Particle  Injection 

The  initial  location  of  each  particle  is  selected  at  random  from  a distribution  function  which 
is  normal  in  the  ^-direction  with  mean  value  equal  to  the  location  of  the  center  of  the  heat 
source  and  variance  equal  to  the  mean  width  of  the  heat  source.  In  the  ^-direction  the 
distribution  function  is  either  normal  with  characteristics  similar  to  that  in  the  ^-direction 
or  exponential.  In  all  cases,  the  distribution  functions  mimic  the  spatial  distribution  of  the 
heat  source.  Also,  particles  are  injected  at  a rate  that  follows  the  rate  of  heat  addition. 

A. 2 Interpolation  of  the  Fields  at  Particles  Locations 

Take  the  2-D  case  and  let  and  77”  be  the  coodinates  of  particle  j,  1 < j < Nv  at  time  tn. 
Let 
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where  [...]  is  the  integer  part  of  the  quantity  in  square  brackets. 

Define  the  function  F as  follows: 

F(Ai,i5  Ai  52  ? ^2,i?  ^2,2  ■>  r,  s)  = Ai^l  — r)(l  — 6)-j-  A2,ir(l  — s)  + — r)s  + A.2,2r<s 
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and  let  Uj,Wj,pj  be  the  velocity  components  and  the  density  respectively  at  the  location  of 
the  jth  particle.  Then  by  bilinear  interpolation: 
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A. 3 Time  Advancement  Scheme  for  Particles 


Now,  each  particle  is  advanced  through  the  differential  equation: 


A ( 5 Vj) 

dt 

J 

drjj 

WjitjiVj) 

dt 

J 

This  corresponds  to  the  differential  equation  describing  the  movement  of  the  particles  in  the 
physical  domain: 


dxj 

dt 

dzj 


dt 


uj{Xj , Zj) 

wj(Xj,Zj) 


The  integration  is  carried  out  in  the  computational  domain  and  the  position  of  the  particle 
is  then  mapped  into  the  physical  domain  for  graphical  purposes. 

The  particles  are  advanced  in  time  according  to  a second  order  Runge-Kutta  time  inte- 
gration scheme  as  follows. 


Then,  with 
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The  new  extrapolated  velocity  components  at  location  -f  htJ'7!?  + are 
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Finally,  the  location  of  the  particle  at  time  level  n + 1 is 


u. 
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The  accuracy  of  this  particle-tracking  routine  was  tested  using  a steady-state  flow  with 
vorticity  in  a rectangular  enclosure.  This  flow  field  was  originally  derived  to  test  the  al- 
gorithms that  integrate  the  fluid  equations  and  is  described  in  [13].  In  the  flow,  particles 
follow  streamlines  which  close  upon  themselves,  a sensitive  test  of  the  quality  of  the  particle- 
tracking integration  scheme. 


B Appendix:  Stability 

The  linear  stability  of  the  lagged-dissipation  method  has  been  checked  for  a single  convection- 
diffusion  equation.  Derivation  of  the  dispersion  relation  is  straightforward.  We  start  with  a 
single  convection- diffusion  equation 

d0_  l^d0_  W dO  _ e 2 

dt  + y/j  d£  + y/J  dr)  J 

where  0 is  the  single  dependent  variable,  U and  V are  the  local  linearized  vector  velocity 
components,  e is  the  dissipation  coefficient,  V2  is  the  Laplacian  and  J is  the  Jacobian  of  the 
spatial  transformation.  We  define 


where  £,■  r)k  = k8r):tn  = n8t. 

Consider  the  difference  form  of  the  convection-diffusion  equation: 
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The  equation  for  the  amplitude  An  becomes 
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Substitution  of  An  — 0n  and  division  by  9n  1 yields  the  dispersion  relation  for  the  lagged- 
diffusion  scheme 

92  A6iT-  [1  -2d(l  -P)]  = 0 

The  roots  of  this  quadratic  equation,  the  amplification  factors  for  the  scheme,  are 


e = 


-iT  ± J-T 2 + 4[1  - 2d(l  - P)] 


If  the  Courant  condition,  < 1,  is  satisfied  and  d = 2cSt}^r  ^ < 1,  then,  for  both  roots, 

\e\  = yji-2d{i  - p ) < i 
since  |P|  < 1,  and  the  scheme  is  stable. 


References 

[1]  Jameson,  A.,  “Computational  Aerodynamics  for  Aircraft  Design,”  Science,  245,  pp. 
361-371,  July  1989. 

[2]  Rehm,  R.G.,  H.C.  Tang,  H.R.  Baum,  J.S.  Sims  and  D.M.  Corley,  “A  Boussinesq  Al- 
gorithm for  Enclosed  Buoyant  Convection  in  Two  Dimensions”,  NIST  Internal  Report, 
NISTIR  4540 , March,  1991. 

[3]  Rehm,  R.G.  and  H.R.  Baum,  “The  Equations  of  Motion  for  Thermally  Driven,  Buoyant 
Flows”,  Journal  of  Research  of  the  NBS,  83,  pp.  297-308,  May- June  1978. 


16 


[4]  Baum,  H.R.,  R.G.  Rehm,  P.D.  Barnett  and  D.M.  Corley,  “Finite  Difference  Calculations 
of  Buoyant  Convection  in  an  Enclosure”,  SIAM  J.  Set.  Stat.  Comput.,  4,  PP-  117-135, 
March  1983. 

[5]  Baum,  H.R.  and  Rehm,  R.G.,  “ Finite  Difference  Solutions  for  Internal  Waves  in  En- 
closures”, SIAM  J.  Sci.  Stat.  Comput.,  5,  No.  4,  pp.  958-977,  1984. 

[6]  Rehm,  R.G.,  P.D.  Barnett,  H.R.  Baum  and  D.M.  Corley,  “Finite  Difference  Calculations 
Of  Buoyant  Convection  in  an  Enclosure:  Verification  of  the  Nonlinear  Algorithm”, 
Applied  Numerical  Mathematics,  1 , pp.  515-529,  North-Holland,  1985. 

[7]  Baum,  H.R.  and  R.G.  Rehm,  “Calculations  of  Three-Dimensional  Buoyant  Plumes  in 
Enclosures”,  Combustion  Science  and  Technology,  40 , pp.  55-77,  Gordon  and  Breach 
Science  Publishers,  1984. 

[8]  Gresho,  P.M.,  “A  Modified  Finite-Element  Method  for  Solving  the  Incompressible 
Navier-Stokes  Equations,”  Large-Scale  Computations  in  Fluid  Mechanics,  AMS  Lec- 
tures in  Applied  Mathematics,  Vol.  22,  Part  1,  pp.  193-240,  1985. 

[9]  Trefethen,  L.N.,  “Numerical  Computation  of  the  Schwarz-Christoffel  Transformation”, 
SIAM  J.  Sci.  Stat.  Comput.,  1 , pp.  82-102,  1980. 

[10]  Simcox,  S.,  Wilkes,  N.S.  and  Jones,  I.P.,  “Fire  at  King’s  Cross  Underground  Station, 
18th  November  1987:  Numerical  Simulation  of  the  Buoyant  Flow  and  Heat  Transfer”, 
Harwell  Report  AERE-G  4677,  May  1988. 

[11]  Cox,  G.,  Chitty,  R.  and  Kumar,  S.,  “Fire  Modeling  and  the  King’s  Cross  Fire  Investi- 
gation”, Letter  to  the  Editor,  Fire  Safety  Journal,  15  pp.  103-106,  1989. 

[12]  Rehm,  R.G.,  Baum,  H.R.,  Lozier,  D.W.,  Tang,  H.C.  and  Sims,  J.,  “Buoyant  Convection 
in  an  Inclined  Enclosure,”  Third  International  Symposium  on  Fire  Safety  Science,  The 
University  of  Edinburgh,  Scotland,  U.K.,  July,  1991. 

[13]  Rehm,  R.G.,  Baum,  H.R.  and  Barnett,  P.D.,  “Buoyant  Convection  Computed  in  a 
Vorticity,  Stream-Function  Formulation,”  Journal  of  Research  of  the  NBS,  87 , No.  2, 
March-April  1982,  pp.  165-185. 


17 


Figure  1.  Window  grid  geometry  demonstrating  the  conformal  transformation  of  the  com- 
putational into  the  physical  domain. 
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Figure  2.  Simulation  of  fire  in  a stair  well,  demonstrating  the  “trench  effect”.  Adiabatic, 
free-slip  boundary  conditions  have  been  imposed.  The  Reynolds  number  is  4 x 104;  the  grid 
size  is  1024  x 256. 
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Figure  3.  Simulation  of  smoke  pouring  from  a window.  Adiabatic,  no-slip  boundary  con- 
ditions have  been  imposed.  The  Reynolds  number  is  1 x 104;  the  grid  size  is  256  x 128. 
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Figure  4.  Comparison  of  four  different  boundary  conditions  for  the  room-window  fire 
scenario.  The  no-slip  condition  represents  the  more  realistic  model,  but  is  limited  by  low 
Reynolds  number. 
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TIME  = 10  0 


Figure  5.  Flow  of  smoke  down  a corridor  with  soffits.  Adiabatic,  no-slip  boundary  condi- 
tions have  been  imposed.  The  Reynolds  number  is  1 x 1 04 ; the  grid  size  is  512  x 128.  This 
geometry  points  out  the  difficulties  faced  by  the  current  methodology  in  handling  the  flow 
around  sharp  corners. 
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Figure  6.  Density  perturbation  in  the  boundary  cells  at  times  t = 16  and  t — 18.  This 
instability  of  the  method  is  not  predicted  by  the  linear  stability  analysis  of  Appendix  B. 
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